notebook.community

Edit and run

See the RSS notebook for introduction.





In [1]:



    


#!wget http://www.remss.com/data/msu/data/netcdf/uat4_tb_v03r03_avrg_chTLT_197812_201308.nc3.nc
#!mv uat4_tb_v03r03_avrg_chTLT_197812_201308.nc3.nc data/



    











In [2]:



    


#!wget http://www.remss.com/data/msu/data/netcdf/uat4_tb_v03r03_avrg_chTMT_197812_201308.nc3.nc
#!mv uat4_tb_v03r03_avrg_chTMT_197812_201308.nc3.nc data/



    











In [3]:



    


#!wget http://www.remss.com/data/msu/data/netcdf/uat4_tb_v03r03_anom_chTLT_197812_201308.nc3.nc
#!mv uat4_tb_v03r03_anom_chTLT_197812_201308.nc3.nc data/

Weight functions





In [4]:



    


%pylab inline
import urllib2
import os
from IPython.display import Image

def download(url, dir):
    """Saves file 'url' into 'dir', unless it already exists."""
    filename = os.path.basename(url)
    fullpath = os.path.join(dir, filename)
    if os.path.exists(fullpath):
        print "Already downloaded:", filename
    else:
        print "Downloading:", filename
        open(fullpath, "w").write(urllib2.urlopen(url).read())



    


















    






Populating the interactive namespace from numpy and matplotlib

















In [5]:



    


download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_TTS.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_TLS.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_tlt_land.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_tlt_ocean.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_tmt_land.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_tmt_ocean.txt", "data")



    


















    






Already downloaded: std_atmosphere_wt_function_chan_TTS.txt
Already downloaded: std_atmosphere_wt_function_chan_TLS.txt
Already downloaded: std_atmosphere_wt_function_chan_tlt_land.txt
Already downloaded: std_atmosphere_wt_function_chan_tlt_ocean.txt
Already downloaded: std_atmosphere_wt_function_chan_tmt_land.txt
Already downloaded: std_atmosphere_wt_function_chan_tmt_ocean.txt

















In [6]:



    


D = loadtxt("data/std_atmosphere_wt_function_chan_TTS.txt", skiprows=6)
h = D[:, 1]
wTTS = D[:, 5]



    











In [7]:



    


D = loadtxt("data/std_atmosphere_wt_function_chan_TLS.txt", skiprows=6)
assert max(abs(h-D[:, 1])) < 1e-12
wTLS = D[:, 5]



    











In [8]:



    


D = loadtxt("data/std_atmosphere_wt_function_chan_tlt_land.txt", skiprows=7)
assert max(abs(h-D[:, 1])) < 1e-12
wTLT_land = D[:, 5]



    











In [9]:



    


D = loadtxt("data/std_atmosphere_wt_function_chan_tlt_ocean.txt", skiprows=7)
assert max(abs(h-D[:, 1])) < 1e-12
wTLT_ocean = D[:, 5]



    











In [10]:



    


D = loadtxt("data/std_atmosphere_wt_function_chan_tmt_land.txt", skiprows=7)
assert max(abs(h-D[:, 1])) < 1e-12
wTMT_land = D[:, 5]



    











In [11]:



    


D = loadtxt("data/std_atmosphere_wt_function_chan_tmt_ocean.txt", skiprows=7)
assert max(abs(h-D[:, 1])) < 1e-12
wTMT_ocean = D[:, 5]



    











In [12]:



    


figure(figsize=(3, 8))
plot(wTLS, h/1000, label="TLS")
plot(wTTS, h/1000, label="TTS")
plot(wTMT_ocean, h/1000, label="TMT ocean")
plot(wTMT_land, h/1000, label="TMT land")
plot(wTLT_ocean, h/1000, label="TLT ocean")
plot(wTLT_land, h/1000, label="TLT land")
xlim([0, 0.2])
ylim([0, 50])
legend()
ylabel("Height [km]")
show()
Image(url="http://www.ssmi.com/msu/img/wt_func_plot_for_web_2012.all_channels2.png", embed=True)



    


















    

















    
Out[12]:








<IPython.core.display.Image at 0x26b6e10>

Netcdf data





In [13]:



    


from netCDF4 import Dataset
from numpy.ma import average



    











In [14]:



    


rootgrp = Dataset('data/uat4_tb_v03r03_avrg_chTLT_197812_201308.nc3.nc')
list(rootgrp.variables)



    


















    
Out[14]:








[u'longitude',
 u'latitude',
 u'time',
 u'climatology_time',
 u'longitude_bounds',
 u'latitude_bounds',
 u'time_bounds',
 u'climatology_time_bounds',
 u'brightness_temperature',
 u'brightness_temperature_climatology',
 u'offset_values',
 u'target_factor_values',
 u'tb_factor_values',
 u'msu_amsu_offsets',
 u'satellites_used']

















In [15]:



    


# 144 values, interval [-180, 180]
longitude = rootgrp.variables["longitude"][:]
# 72 values, interval [-90, 90]
latitude = rootgrp.variables["latitude"][:]
# 144 rows of [min, max]
longitude_bounds = rootgrp.variables["longitude_bounds"][:]
# 72 rows of [min, max]
latitude_bounds = rootgrp.variables["latitude_bounds"][:]
# time in days, 1978 - today
time = rootgrp.variables["time"][:]
# time in years
years = time / 365.242 + 1978
# 12 values: time in days for 12 months in a year
time_climatology = rootgrp.variables["climatology_time"][:]
# (time, latitude, longitude)
brightness_temperature = rootgrp.variables["brightness_temperature"][:]
# (time_climatology, latitude, longitude)
brightness_temperature_climatology = rootgrp.variables["brightness_temperature_climatology"][:]

The weighting function (based on area) is given by: $$ w_\theta = \sin{\pi\over144}\cos\theta $$ $$ \sum_\theta w_\theta = 1 $$





In [16]:



    


w_theta = sin(pi/144) * cos(latitude*pi/180)
sum(w_theta)



    


















    
Out[16]:








0.99999988

















In [17]:



    


Tavg = average(brightness_temperature, axis=2)
Tavg = average(Tavg, axis=1, weights=w_theta)
plot(years, Tavg-273.15)
xlabel("Year")
ylabel("T [C]")
title("TLT (Temperature Lower Troposphere)")
show()

The temperature oscillates each year. To calculate the "anomaly", we subtract from each month its average temperature:





In [18]:



    


Tanom = empty(Tavg.shape)
for i in range(12):
    Tanom[i::12] = Tavg[i::12] - average(Tavg[i::12])

We calculate linear fit





In [19]:



    


from scipy.optimize import curve_fit
# Skip the first year, start from 1979, that's why you see the "12" here and below:
Y0 = years[12]
def f(x, a, b):
    return a*(x-Y0)+b
(a, b), pcov = curve_fit(f, years[12:], Tanom[12:])
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par              dev"
print a, adev
print b, bdev



    


















    






par              dev
0.0128077506785 0.000847402784383
-0.220550651996 0.0169189854372

And compare against official graph + trend. As can be seen, the agreement is perfect:





In [20]:



    


from matplotlib.ticker import MultipleLocator
figure(figsize=(6.6, 3.5))
plot(years, Tanom, "b-", lw=0.7)
plot(years, a*(years-Y0)+b, "b-", lw=0.7, label="Trend = $%.3f \pm %.3f$ K/decade" % (a*10, adev*10))
xlim([1979, 2014])
ylim([-1.2, 1.2])
gca().xaxis.set_minor_locator(MultipleLocator(1))
legend()
xlabel("Year")
ylabel("Temperature Anomaly [K]")
title("TLT (Temperature Lower Troposphere)")
show()
Image(url="http://www.remss.com/data/msu/graphics/TLT/plots/RSS_TS_channel_TLT_Global_Land_And_Sea_v03_3.png", embed=True)



    


















    

















    
Out[20]:








<IPython.core.display.Image at 0x30f3dd0>

Maps

Let's show some plots.





In [21]:



    


from mpl_toolkits.basemap import Basemap
X, Y = meshgrid(longitude, latitude)



    











In [22]:



    


idx = 12*30
#X, Y = meshgrid(concatenate([longitude, [longitude[0]]]), concatenate([latitude, [latitude[0]]]))
N = 12

print 'Temperatures (%s)' % int(years[idx])
Tall=(brightness_temperature[idx:idx+N, :, :]-273.15)
vmin = Tall.min()
vmax = Tall.max()

from mpl_toolkits.basemap import Basemap
for i in range(N):
    T = brightness_temperature[idx+i, :, :]-273.15

    figure(figsize=(18, 6))
    
    subplot(1, 3, 1)
    m = Basemap(projection='ortho',lat_0=45, lon_0=-100, resolution='l')
    m.drawcoastlines()
    m.pcolor(X, Y, T, latlon=True, vmin=vmin, vmax=vmax)
    
    subplot(1, 3, 2)
    m = Basemap(projection='hammer',lon_0=0)
    m.drawcoastlines()
    p = m.pcolor(X, Y, T, latlon=True, vmin=vmin, vmax=vmax)
    cbar = m.colorbar(p, location='bottom', pad="5%")
    cbar.set_label('$^\circ C$')
    
    subplot(1, 3, 3)
    m = Basemap(llcrnrlon=0,llcrnrlat=-90,urcrnrlon=360,urcrnrlat=90,projection='mill')
    m.drawcoastlines()
    m.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
    m.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1])
    m.pcolor(X, Y, T, latlon=True, vmin=vmin, vmax=vmax)

    suptitle("Month: %d" % (i+1))

    show()



    


















    






Temperatures (2008)

Remove seasonality (=calculate anomaly):





In [23]:



    


T = brightness_temperature
T2 = empty(T.shape)
for i in range(12):
    T2[i::12, :, :] = T[i::12, :, :] - average(T[i::12, :, :], axis=0)



    











In [24]:



    


m = Basemap(projection='hammer',lon_0=0)
m.drawcoastlines()
p = m.pcolor(X, Y, T2[12*20, :, :], latlon=True)
cbar = m.colorbar(p, location='bottom', pad="5%")
cbar.set_label('$^\circ C$')

Calculate trends:





In [25]:



    


import numpy.ma as ma
def trend(years, series):
    Y0 = years[12]
    def f(x, a, b):
        return a*(x-Y0)+b
    (a, b), pcov = curve_fit(f, years[12:], series[12:])
    if isinstance(pcov, float):
        return ma.masked
    adev = sqrt(pcov[0, 0])
    bdev = sqrt(pcov[1, 1])
    return a*10



    











In [26]:



    


T3 = ma.empty(T2.shape[1:])
print T3.shape
for i in range(size(T2, 1)):
    print i
    for j in range(size(T2, 2)):
        T3[i, j] = trend(years, T2[:, i, j])



    


















    






(72, 144)
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

















In [27]:



    


m = Basemap(projection='hammer',lon_0=0)
m.drawcoastlines()
p = m.pcolor(X, Y, T3, latlon=True, vmin=-0.6, vmax=0.6)
cbar = m.colorbar(p, location='bottom', pad="5%")
cbar.set_label('Trend [$^\circ$C/decade]')

Average over longitude:





In [28]:



    


T4_TLT = average(T3, axis=1)



    











In [29]:



    


plot(latitude, T4_TLT)



    


















    
Out[29]:








[<matplotlib.lines.Line2D at 0x52a6450>]










    
























In [30]:



    


rootgrp = Dataset('data/uat4_tb_v03r03_avrg_chTMT_197812_201308.nc3.nc')
brightness_temperature = rootgrp.variables["brightness_temperature"][:]
T = brightness_temperature
T2 = empty(T.shape)
for i in range(12):
    T2[i::12, :, :] = T[i::12, :, :] - average(T[i::12, :, :], axis=0)
T3 = ma.empty(T2.shape[1:])
print T3.shape
for i in range(size(T2, 1)):
    print i
    for j in range(size(T2, 2)):
        T3[i, j] = trend(years, T2[:, i, j])
T4_TMT = average(T3, axis=1)



    


















    






(72, 144)
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

















In [31]:



    


rootgrp = Dataset('data/uat4_tb_v03r03_avrg_chTLS_197812_201308.nc3.nc')
brightness_temperature = rootgrp.variables["brightness_temperature"][:]
T = brightness_temperature
T2 = empty(T.shape)
for i in range(12):
    T2[i::12, :, :] = T[i::12, :, :] - average(T[i::12, :, :], axis=0)
T3 = ma.empty(T2.shape[1:])
print T3.shape
for i in range(size(T2, 1)):
    print i
    for j in range(size(T2, 2)):
        T3[i, j] = trend(years, T2[:, i, j])
T4_TLS = average(T3, axis=1)



    


















    






(72, 144)
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

















In [32]:



    


profile= wTLT_ocean/max(wTLT_ocean) * T4_TLT[:, None]+wTMT_ocean/max(wTMT_ocean) * T4_TMT[:, None]+wTLS/max(wTLS) * T4_TLS[:, None]

This should be compared with Fig. 2H in:

Santer, B. D., Painter, J. F., Bonfils, C., Mears, C. a., Solomon, S., Wigley, T. M. L., … Wentz, F. J. (2013). Human and natural influences on the changing thermal structure of the atmosphere. Proceedings of the National Academy of Sciences, 6–11. doi:10.1073/pnas.1305332110





In [33]:



    


imshow(profile[:, :50], vmin=-0.3, vmax=0.3)
colorbar()



    


















    
Out[33]:








<matplotlib.colorbar.Colorbar instance at 0x3406e60>

Content source: certik/climate

Similar notebooks:

notebook.community | gallery | about